#Setup
library(Seurat)
library(tidyverse)
library(scales)
library(viridis)
library(patchwork)
library(ggrepel)
#set working directory to the location of the file
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
source("../my_functions.R") #Loading a couple of custom functions for plotting
b_cells <- read.table('b_cells.txt') %>% pull(V1)
#-Subset B cell clusters
#Idents(Neurogenic_lineage) = "SCT_snn_res.2"
Head <- subset(integrated_exp, cells = b_cells)
Head <- Head %>%
SCTransform() %>%
RunPCA(assay = "SCT", npcs = 100) %>%
IntegrateLayers(method = CCAIntegration,
normalization.method = "SCT",
verbose = F) %>%
RunUMAP(dims = 1:100, reduction = "integrated.dr") %>%
FindNeighbors(reduction = "integrated.dr", dims = 1:100) %>%
FindClusters(resolution = seq(0.5, 2, 0.5), graph.name = "SCT_snn")
DimPlot(Head, label=T, shuffle = T, group.by ="SCT_snn_res.1") +
NoLegend() +
coord_fixed() +
NoAxes()
umap_original <- Head@reductions$umap@cell.embeddings
# Custom function to create a 2D rotation matrix
rot2 <- function(angle_radians) {
matrix(c(cos(angle_radians), -sin(angle_radians),
sin(angle_radians), cos(angle_radians)), nrow = 2)
}
# Function to rotate multiple points using a matrix of coordinates
rotate_points <- function(points, angle_degrees) {
# Convert the angle from degrees to radians
angle_radians <- angle_degrees * pi / 180
# Get the rotation matrix
rotation_matrix <- rot2(angle_radians)
# Apply the rotation matrix to the matrix of points (transpose for correct multiplication)
rotated_points <- t(rotation_matrix %*% t(points))
return(rotated_points)
}
# Rotate all points by 45 degrees
umap_rotated <- rotate_points(umap_original, 90)
colnames(umap_rotated) <- c("umap_1", "umap_2")
Head@reductions$umap@cell.embeddings <- umap_rotated
DimPlot(Head, label=T, shuffle = T, group.by ="SCT_snn_res.1") +
NoLegend() +
coord_fixed() +
NoAxes()
DimPlot(Head, group.by = c("cell_type_label", "region")) & NoAxes() & coord_fixed()
FeaturePlot(Head, "Wpre", order = T) + scale_color_viridis(option = "A")
####-Identifying Tdtomato+ (Wpre+) cells
VlnPlot(Head, "Wpre", group.by = "batch", split.by = "region") #There is one cell with Wpre > 2 in the Wedge region. This cell could be a contaminant or a B1 cell that transitioned to B2 between injection and cell dissociation.
wpre.data %>% filter(batch == "batch_2" & region == 'LW') %>% pull(wpre_expression) %>% quantile(probs = 0.875)
87.5%
1.098612
#Selecting Wpre+ cells, defined as cells with normalized Wpre expression > 1.39, excluding cells coming from the Wedge region.
wpre.data <- wpre.data %>%
mutate(wpre_label = case_when(batch == "batch_2" &
wpre_expression > batch2_thres &
region!= "Wedge" ~ "wpre+",
batch == "batch_3" &
wpre_expression > batch3_thres &
region != "Wedge" ~ "wpre+",
.default = "wpre-"))
#Calculating the number of unlabeled B1 cells in the dataset:
wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch, wpre_label) %>%
summarise(n = n()) %>%
left_join(wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch) %>%
summarize(batch_n = n()), by = "batch") %>%
mutate(pct = n/batch_n)
####-Identifying non-infected B1 cells
#Assuming that these cells are the most similar to Wpre+ cells.
#Splitting our dataset into different batches, selecting only cells that are not in the Wedge region.
wpre_cells_batch2 <- wpre.data %>%
filter(wpre_label == "wpre+" &
batch == "batch_2") %>%
rownames()
wpre_cells_batch3 <- wpre.data %>%
filter(wpre_label == "wpre+" &
batch == "batch_3") %>%
rownames()
nonwedge_batch2_cells <- Head@meta.data %>% filter(region != "Wedge" & batch == "batch_2") %>% rownames()
nonwedge_batch3_cells <- Head@meta.data %>% filter(region != "Wedge" & batch == "batch_3") %>% rownames()
Head_batch2 <- Head %>% subset(cells = nonwedge_batch2_cells) %>%
SCTransform() %>%
RunPCA(assay = "SCT",
npcs = 100) %>%
FindNeighbors(dims = 1:100,
k.param = 30,
return.neighbor = T)
Head_batch3 <- Head %>% subset(cells = nonwedge_batch3_cells) %>%
SCTransform() %>%
RunPCA(assay = "SCT",
npcs = 100) %>%
FindNeighbors(dims = 1:100,
k.param = 30,
return.neighbor = T)
#Inferring unlabeled b1 cells in batch 2:
wpre_index_batch2 <- which(Head_batch2@neighbors$SCT.nn@cell.names %in% wpre_cells_batch2)
good.neighbors_batch2 <- c()
bad.neighbors <- wpre_index_batch2
for (i in wpre_index_batch2) {
all.neighbors <- Head_batch2@neighbors$SCT.nn@nn.idx[i,2:30]
good.neighbors_batch2 <- c(good.neighbors_batch2,
(all.neighbors[!all.neighbors %in% bad.neighbors])[1:3])
bad.neighbors <- c(bad.neighbors, good.neighbors_batch2) %>% unique()
}
wpre_neighbors_batch2 <- Head_batch2@neighbors$SCT.nn@cell.names[good.neighbors_batch2]
wpre_neighbors_batch2 <- wpre_neighbors_batch2[!is.na(wpre_neighbors_batch2)]
#Inferring unlabeled b1 cells in batch 3:
wpre_index_batch3 <- which(Head_batch3@neighbors$SCT.nn@cell.names %in% wpre_cells_batch3)
good.neighbors_batch3 <- c()
bad.neighbors <- wpre_index_batch3
for (i in wpre_index_batch3) {
all.neighbors <- Head_batch3@neighbors$SCT.nn@nn.idx[i,2:30]
good.neighbors_batch3 <- c(good.neighbors_batch3, (setdiff(all.neighbors, bad.neighbors))[1:3])
bad.neighbors <- c(bad.neighbors, good.neighbors_batch3)
}
wpre_neighbors_batch3 <- Head_batch3@neighbors$SCT.nn@cell.names[good.neighbors_batch3]
wpre_neighbors_batch3 <- wpre_neighbors_batch3[!is.na(wpre_neighbors_batch3)]
wpre.data[, "b_type"] <- "b2"
wpre.data[c(wpre_neighbors_batch2, wpre_neighbors_batch3),"b_type"] <- "b1 - nn"
wpre.data[c(wpre_cells_batch2, wpre_cells_batch3), "b_type"] <- "b1 - labeled"
wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch, b_type) %>%
summarise(n = n()) %>%
left_join(wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch) %>%
summarize(batch_n = n()), by = "batch") %>%
mutate(pct = n/batch_n)
Head@meta.data[, "b_type"] <- "b2"
Head@meta.data[c(wpre_neighbors_batch2, wpre_neighbors_batch3),"b_type"] <- "b1 - nn"
Head@meta.data[c(wpre_cells_batch2, wpre_cells_batch3), "b_type"] <- "b1 - labeled"
DimPlot(Head, group.by = "b_type", shuffle = T) & NoAxes() & coord_fixed()
DimPlot(Head, group.by = "b_type", split.by = "b_type") & NoAxes() & coord_fixed()
####-Label transfer to uninfected cells
Head@active.assay = "RNA"
non_infected_features <- rownames(Head)[!rownames(Head) %in% c('Wpre', 'GFP', 'Cre')]
Head_batch1 <- subset(Head, subset = batch == "batch_1", features = non_infected_features) %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA(assay = "RNA",
npcs = 20) %>%
FindNeighbors(dims = 1:20)
Head_batch23 <- subset(Head, subset = batch %in% c("batch_2", "batch_3"))
Head_batch23 <- Head_batch23 %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA(assay = "RNA",
npcs = 20)
Head_batch23 <- IntegrateLayers(Head_batch23, method = CCAIntegration, orig.reduction = "pca", dims = 1:20)
Head.anchors <- Head_batch23 %>%
FindTransferAnchors(reference = .,
query = Head_batch1,
dims = 1:20,
reference.reduction = "pca")
predictions <- TransferData(anchorset = Head.anchors, refdata = Head_batch23$b_type, dims = 1:20)
inferredb1_batch1_cells <- predictions %>% filter(predicted.id != "b2") %>% rownames()
Head@meta.data[inferredb1_batch1_cells, "b_type"] <- "b1 - nn"
Head@meta.data %>%
filter(region != "Wedge") %>%
group_by(batch, b_type) %>%
summarise(n = n()) %>%
left_join(wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch) %>%
summarize(batch_n = n()), by = "batch") %>%
mutate(pct = n/batch_n)
FeaturePlot(Head, c("S100a6", "Egfr", "Crym", "Urah", "Tfap2c"), order = T, ncol = 5) &
scale_color_viridis(option = "A") &
coord_fixed() & simple
Idents(Head) = "SCT_snn_res.1"
Head@meta.data <- Head@meta.data %>%
mutate(cell_subtype = case_when(
SCT_snn_res.1 %in% c(4, 2, 8, 11, 10) ~ "aB cells",
SCT_snn_res.1 == 0 ~ "Ventral Subpallial B cells",
SCT_snn_res.1 %in% c(3, 5) ~ "Dorsal Pallial B cells",
SCT_snn_res.1 %in% c(6, 1, 9, 7) ~ "Dorsal Subpallial B cells",
TRUE ~ "Unknown"
))
Head@meta.data <- Head@meta.data %>%
mutate(activation_state = case_when(
cell_subtype == "aB cells" ~ "activated",
TRUE ~ "quiescent"
))
Idents(Head) = "SCT_snn_res.1"
DimPlot(Head,
group.by = "cell_subtype",
label=F,
cols = c("cyan","blue","dodgerblue2","dodgerblue4", "dodgerblue3")) +
NoLegend() +
coord_fixed() +
NoAxes()
DimPlot(Head, group.by ="region", label = F, shuffle = F, cols = c("grey", "purple", "dodgerblue2"), pt.size = 1.5) & NoAxes() & coord_fixed()
#B1 score -----------------------------
Head@active.assay <- "SCT"
b.markers.list = list(b1 = c("Atf3", "Riiad1", "Foxj1", "Gadd45b", "Tagln2", "Emp1"))
Head = AddModuleScore(Head, features = b.markers.list)
Head@meta.data = Head@meta.data %>% dplyr::rename(B1_score = Cluster1)
#Re-scaling Module Scores:These lines rescale the module scores (B1_score and B2_score) so that they have a mean of 0 and a standard deviation of 1. This step standardizes the scores, making them comparable across different modules.
Head$B1_score = Head$B1_score %>% rescale()
FeaturePlot(Head, c("B1_score"), order = T, pt.size = 1.5) &
scale_color_viridis(option = "magma") &
NoAxes() &
coord_fixed()
saveRDS(Head, "data/Head.rds")
integrated_exp <- readRDS("../pre-processing/integrated_exp_withp365.rds")
mDimPlot(integrated_exp, group.by = "cell_type", label = T, shuffle = T, legend = F, repel = T)
b_cells <- integrated_exp@meta.data %>% filter(cell_type == 'B cells') %>% rownames()
Head <- subset(integrated_exp, cells = b_cells)
Head <- Head %>%
SCTransform() %>%
RunPCA(assay = "SCT", npcs = 100) %>%
IntegrateLayers(method = CCAIntegration,
normalization.method = "SCT",
verbose = F) %>%
RunUMAP(dims = 1:100, reduction = "integrated.dr") %>%
FindNeighbors(reduction = "integrated.dr", dims = 1:100) %>%
FindClusters(resolution = seq(0.5, 2, 0.5), graph.name = "SCT_snn")
Running SCTransform on assay: RNA
Running SCTransform on layer: counts.1
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 18667 by 3092
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 3092 cells
Found 67 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 18667 genes
Computing corrected count matrix for 18667 genes
Calculating gene attributes
Wall clock passed: Time difference of 40.20527 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|==================== | 25%
|
|======================================== | 50%
|
|============================================================ | 75%
|
|================================================================================| 100%
Running SCTransform on layer: counts.2
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 16588 by 520
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 520 cells
Found 19 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 16588 genes
Computing corrected count matrix for 16588 genes
Calculating gene attributes
Wall clock passed: Time difference of 10.79082 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|==================== | 25%
|
|======================================== | 50%
|
|============================================================ | 75%
|
|================================================================================| 100%
Running SCTransform on layer: counts.3
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 17437 by 2293
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 2293 cells
Found 28 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 17437 genes
Computing corrected count matrix for 17437 genes
Calculating gene attributes
Wall clock passed: Time difference of 25.02378 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|==================== | 25%
|
|======================================== | 50%
|
|============================================================ | 75%
|
|================================================================================| 100%
Centering data matrix
|
| | 0%
|
|=========== | 14%
|
|======================= | 29%
|
|================================== | 43%
|
|============================================== | 57%
|
|========================================================= | 71%
|
|===================================================================== | 86%
|
|================================================================================| 100%
Centering data matrix
|
| | 0%
|
|=========== | 14%
|
|======================= | 29%
|
|================================== | 43%
|
|============================================== | 57%
|
|========================================================= | 71%
|
|===================================================================== | 86%
|
|================================================================================| 100%
Centering data matrix
|
| | 0%
|
|=========== | 14%
|
|======================= | 29%
|
|================================== | 43%
|
|============================================== | 57%
|
|========================================================= | 71%
|
|===================================================================== | 86%
|
|================================================================================| 100%
Warning: Different cells and/or features from existing assay SCTSet default assay to SCT
Warning: The following 64 features requested have not been scaled (running reduction without them): Ifi27l2a, GFP, Cre, Lcn2, C3, Wpre, Ifi44, Gm42047, H2-Aa, Phf11d, Cfb, H2-Eb1, Slfn8, Tgm1, Gm11728, AW112010, BC028528, Olfr1369-ps1, Fam216b, Gm10457, 4930477G07Rik, Cxcl5, Ssc4d, Shoc1, Oasl1, Oas1a, Tex21, Tgtp2, Hs3st4, Mlkl, Wfdc2, Apol9a, Cxcl9, Zyg11a, Myoc, C2cd4b, Gm10636, Tbx20, Actn3, Gm13538, Gm13219, Cdcp1, Ccl5, Btbd16, Acod1, Npffr1, Scgb1c1, Gm16090, Atp12a, Syce3, Trarg1, Gm29773, Ccl19, Sfta2, Ifi204, 4930445E18Rik, Gm32913, Slc28a1, Tbc1d30, 4732419C18Rik, Gm12185, Pde6c, Fgf23, Slc44a4PC_ 1
Positive: Egfr, Nrxn3, Kcnh7, Gm29260, Tmsb4x, Slit2, Ccnd2, Fabp7, Ptprz1, Lrrc7
Pak3, Nol4, Sox11, Nrxn1, Dscaml1, Zbtb20, Setbp1, Sv2c, Sox2ot, Tcf12
Slc4a7, Mir99ahg, Zeb1, Msi2, Nrg1, Marcks, Map2, Sox4, Auts2, Shc3
Negative: Mt1, Junb, Klf2, Clu, Crym, Cebpd, Mt3, Gfap, S100a6, Id3
Gadd45g, Zfp36, Atf3, Ccn1, Fth1, Cebpb, Cst3, Cxcl10, Ier2, Fxyd1
Mt2, Btg2, Apoe, Aldoc, Prdx6, Dusp1, Klf4, Ifitm3, Fos, Rgcc
PC_ 2
Positive: Crym, Klf2, Cxcl10, Atf3, Dnajb1, Rsad2, Cebpd, Egr3, Isg15, Emp1
Klf6, Jun, Egr2, Klf4, Gbp2, Fas, Zfp36, Serpina3n, Nfatc2, Arc
H2-K1, Ddit3, Actb, Ncam2, Junb, Prss56, Ier2, Homer1, Btg2, Slit2
Negative: Nrg1, Pla2g7, Sparcl1, Mt1, Clu, Aldoc, Thbs4, Ntrk2, Gpc5, Gfap
Igfbp5, Gm29260, Gabrb1, Dcc, Mdga2, Slc15a2, Pax6, Fxyd6, Cst3, Rmst
Cdk11b, Cldn10, Slc1a2, Lama2, Lrp4, Acsl3, Cadm2, Apoe, Fbxo2, Gm42418
PC_ 3
Positive: Hes5, Btg2, Homer1, Jun, Fos, Ier2, Cdk11b, Pde10a, Jund, Junb
Hes1, Trpm3, Dnajb1, Klf4, Egr1, Arih1, Nr4a3, Nrg1, Nr4a1, Dusp1
Vps37b, Ddit3, Fosb, Gpc5, Urah, Usp2, Coq10b, Gm39325, Nkain2, Maff
Negative: Mt1, Gfap, Basp1, Nnat, Tmsb10, Sparcl1, Dlx6os1, Sox2ot, Stmn2, Nav3
Aqp4, S100a6, Nrxn3, Clu, Tubb3, Ifitm3, Mt3, Mt2, Stmn1, Dcx
Gm42418, Fth1, Celf4, Tuba1a, Cebpd, Igfbpl1, Lrrc7, Serpina3n, Tubb5, Adarb2
PC_ 4
Positive: Crym, Frmd4a, Klf2, Grm3, Nrxn1, Ntm, Ncam2, Dpyd, Ptprt, Gabrb1
Adamts18, Lsamp, Mir99ahg, Prex2, Grid2, Lrrc4c, B3galt1, Slit2, Cadm2, Sgcd
Rora, Adgrl3, Hdac9, Negr1, Pcdh9, Fars2, Cebpd, Prss56, Slc1a2, Zbtb20
Negative: Cdk11b, Pde10a, Gfap, Vps37b, Tmsb4x, Emd, Mt1, Slc25a3, Urah, Arih1
Dad1, Rps24, Tubb5, Nop53, Rpl13, Snhg15, Cebpb, Rpsa, Rplp1, Trmt61b
Ccnd2, Rpl32, Rps8, Maff, Rpl23, Rps27a, Nr4a2, Nrxn3, Rplp0, Nme2
PC_ 5
Positive: Dlx6os1, Nav3, Stmn2, Nrxn3, Sox2ot, Tubb3, Dcx, Basp1, Celf4, Adarb2
Erbb4, Gm38505, Igfbpl1, Gad2, Cacna1c, Tmsb10, Junb, Anks1b, Shtn1, Sp9
Btg2, Enox1, Plcl1, 2610307P16Rik, Myt1l, Map1b, Nol4, Hes5, Lrrc7, Ncam1
Negative: Fabp7, Egfr, Sv2c, Rpsa, Rplp1, Rps8, Rpl13, Nme2, Rplp0, Rpl32
Rpl39, Lrrfip1, Shc3, Rpl23, Rps24, Rps27a, Rps12, Slc4a7, Emp1, Ptprz1
Tmsb4x, Eef1a1, Ptx3, Gm29260, Hsp90ab1, Pvt1, Ankrd28, Tnfrsf12a, Rpl15, Rorb
| | 0 % ~calculating
|+++++++++++++++++ | 33% ~00s
|++++++++++++++++++++++++++++++++++ | 67% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s
| | 0 % ~calculating
|+++++++++++++++++ | 33% ~09s
|++++++++++++++++++++++++++++++++++ | 67% ~05s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=15s
13:06:11 UMAP embedding parameters a = 0.9922 b = 1.112
13:06:11 Read 5905 rows and found 100 numeric columns
13:06:11 Using Annoy for neighbor search, n_neighbors = 30
13:06:11 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:06:12 Writing NN index file to temp file /tmp/RtmpvrQ0a0/file1b8b1f474deaf1
13:06:12 Searching Annoy index using 1 thread, search_k = 3000
13:06:14 Annoy recall = 100%
13:06:19 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
13:06:26 Initializing from normalized Laplacian + noise (using RSpectra)
13:06:26 Commencing optimization for 500 epochs, with 263284 positive edges
Using method 'umap'
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:06:36 Optimization finished
Computing nearest neighbor graph
Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 5905
Number of edges: 415420
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8226
Number of communities: 8
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 5905
Number of edges: 415420
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.7517
Number of communities: 11
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 5905
Number of edges: 415420
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.7033
Number of communities: 15
Elapsed time: 0 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 5905
Number of edges: 415420
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.6588
Number of communities: 16
Elapsed time: 0 seconds
mDimPlot(Head, label=T, shuffle = T, group.by ="SCT_snn_res.1", legend = F)
umap_original <- Head@reductions$umap@cell.embeddings
# Rotate all points by 45 degrees
umap_rotated <- rotate_points(umap_original, 180)
colnames(umap_rotated) <- c("umap_1", "umap_2")
Head@reductions$umap@cell.embeddings <- umap_rotated
mDimPlot(Head, label=T, shuffle = T, group.by ="SCT_snn_res.0.5", legend = F)
#- B1 / B2 identification ### Identifying Tdtomato+ (Wpre+) cells
#Assuming that B1 and B2 cells co-exist in similar numbers at P30, and only 25% of B1 cells are infected, we should select the top 12.5% Wpre-expressing cells in each batch.
wpre.data = Head@meta.data %>%
mutate(wpre_expression = Head@assays$SCT@data["Wpre",] %>%
as.data.frame() %>%
pull("."))
batch2_thres <- wpre.data %>%
filter(batch == "batch_2" &
region == 'LW' &
age == 'p30') %>%
pull(wpre_expression) %>%
quantile(probs = 0.875)
batch3_thres <- wpre.data %>%
filter(batch == "batch_3" &
region == 'LW' &
age == 'p30') %>%
pull(wpre_expression) %>%
quantile(probs = 0.875)
batch2_thres
87.5%
1.386294
batch3_thres
87.5%
1.609438
VlnPlot(Head, "Wpre", group.by = "batch", split.by = "region") + geom_hline(yintercept = c(batch2_thres, batch3_thres, 2.1)) #There is one cell with Wpre > 2 in the Wedge region. This cell could be a contaminant or a B1 cell that transitioned to B2 between injection and cell dissociation.
VlnPlot(Head, "Wpre", group.by = "age", split.by = "region") + geom_hline(yintercept = c(batch2_thres, batch3_thres, 2.1)) #There is one cell with Wpre > 2 in the Wedge region. This cell could be a contaminant or a B1 cell that transitioned to B2 between injection and cell dissociation.
#Selecting Wpre+ cells, defined as cells with high Wpre expression, excluding cells coming from the Wedge region
wpre.data <- wpre.data %>%
mutate(wpre_label = case_when(age != 'p30' | batch == 'batch_1' ~ NA,
batch == "batch_2" &
age == 'p30' &
wpre_expression >= batch2_thres &
region!= "Wedge" &
bcell_subtype != "Dorsal Pallial B cells" ~
"wpre+",
batch == "batch_3" &
age == 'p30' &
wpre_expression >= batch3_thres &
region != "Wedge" &
bcell_subtype != "Dorsal Pallial B cells" ~
"wpre+",
.default = "wpre-"))
#Calculating the number of unlabeled B1 cells in the dataset:
wpre.data %>%
filter(region != "Wedge"
) %>%
group_by(batch, wpre_label) %>%
summarise(n = n()) %>%
left_join(wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch) %>%
summarize(batch_n = n()), by = "batch") %>%
mutate(pct = n/batch_n)
`summarise()` has grouped output by 'batch'. You can override using the `.groups` argument.
#Assuming that these cells are the most similar to Wpre+ cells.
#Splitting our dataset into different batches, selecting only cells that are not in the Wedge region.
wpre_cells_batch2 <- wpre.data %>%
filter(wpre_label == "wpre+" &
batch == "batch_2") %>%
rownames()
wpre_cells_batch3 <- wpre.data %>%
filter(wpre_label == "wpre+" &
batch == "batch_3") %>%
rownames()
batch2_lw_cells <- Head@meta.data %>% filter(region == "LW" &
batch == "batch_2") %>% rownames()
batch3_lw_cells <- Head@meta.data %>% filter(region == "LW" &
batch == "batch_3") %>% rownames()
Head_batch2 <- Head %>% subset(cells = batch2_lw_cells) %>%
SCTransform() %>%
RunPCA(assay = "SCT",
npcs = 100) %>%
FindNeighbors(dims = 1:10,
k.param = 30,
return.neighbor = T)
Running SCTransform on assay: RNA
Running SCTransform on layer: counts.2
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 16203 by 434
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 434 cells
Found 35 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 16203 genes
Computing corrected count matrix for 16203 genes
Calculating gene attributes
Wall clock passed: Time difference of 9.660439 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|==================== | 25%
|
|======================================== | 50%
|
|============================================================ | 75%
|
|================================================================================| 100%
Warning: Different cells and/or features from existing assay SCTSet default assay to SCT
PC_ 1
Positive: Egfr, Kcnh7, Marcks, Tmsb4x, Ptprz1, Pak3, Ccnd2, Magi1, Sall3, Fgd4
Fabp7, Shc3, Dscaml1, Slc4a7, Tcf12, Ppp1r14b, Map2, Hnrnpab, Hmgn2, Sox11
Grin2b, Slit2, H2afz, Msi2, Rpl32, Aff3, Set, Nell2, Gm29260, Nol4
Negative: Clu, Mt3, Gfap, Mt1, S100a6, Fxyd1, Junb, Aldoc, Fth1, Ifitm3
Cebpd, H2-K1, Apoe, Cst3, Ifit3, Zfp36, Serping1, B2m, Sat1, Sparc
Ccn1, Rarres2, Atf3, Id3, Klf2, Ifit3b, H2-D1, Dusp1, Sdc4, Prdx6
PC_ 2
Positive: Nrg1, Ntrk2, Fos, 6330403K07Rik, Acsl3, Pla2g7, Grid2, Rgs6, Gpc5, Mt1
Gabrb1, Btg2, Lsamp, Thbs4, Rmst, Lama2, Cadm2, Dusp1, Slc15a2, Dcc
A330093E20Rik, Slc1a2, Htra1, Sgip1, Kcnn2, Rora, Tsc22d3, Cst3, Usp2, Lrp4
Negative: Bst2, B2m, Ly6e, Ifitm3, Isg15, H2-D1, Iigp1, Irf7, Ifit3, Gbp2
H2-K1, Rsad2, Cxcl10, Ifit3b, Psmb8, Oasl2, H2-T23, Rnf213, Igtp, Ifi44
Ifit1, Gm4951, Serping1, Tspo, Emp1, Psmb10, Irgm1, Vim, Psme1, Crym
PC_ 3
Positive: Nrg1, Gm29260, Sparcl1, Pla2g7, Epha5, Dlgap1, Ntrk2, Rlbp1, Sema6d, Csmd1
Pax6, Kcnk10, Gfap, Clu, Mgat5, Hopx, Rgs6, Limch1, Ptprz1, Pls1
Luzp2, Lrp4, Adam23, Dgki, Mir99ahg, Add3, Mdga2, Aldoc, Arx, Astn2
Negative: Btg2, Junb, Dnajb1, Egr2, Hes5, Fos, Homer1, Klf2, Nr4a1, Jund
Ddit3, Zfp36, Maff, Egr3, Ier2, Klf4, Rhob, Ftl1, Crym, Jun
Notum, Actg1, Fosb, Idi1, Atf3, Egr1, Cebpb, Ubald1, Prss56, Dusp1
PC_ 4
Positive: Dclk1, Nr4a3, Pde10a, Nrg1, Pla2g7, Mt1, Ntrk2, Gfap, Cebpb, Lmna
Clu, Gsn, Slc15a2, Ddx3y, Tuba1b, Srxn1, Emd, Mt2, Usp2, Igfbp5
Slc25a3, Neat1, Gm10260, Coq10b, Rorb, Emp1, Prm1, H2afz, Bag3, S100a10
Negative: Crym, Hes5, Hes1, Frmd4a, Dpyd, Adamts18, Grm3, Rbms3, Lrrc4c, Nrxn1
Prex2, Notum, Gm10561, Adgrl3, Sgcd, Egr2, Ncam1, Prss56, Pbx1, Ptprt
Ntm, Xist, Sntb1, Ncam2, Sema6a, Sox2ot, Negr1, Celf2, Slc1a3, Grik2
PC_ 5
Positive: Ifrd1, Zswim6, Homer1, Emp1, March3, Nfatc2, Sik2, Lmna, Rorb, Creb5
Dclk1, Hmgcr, Cebpb, Zbtb11, Idi1, Fosb, Nebl, Ascc3, Insig1, Negr1
Ddit3, Msmo1, Maff, Neat1, Trib1, Dleu2, Nr4a3, Nr4a1, Atf3, Ankrd28
Negative: Gadd45g, Ckb, Ifi27, Itm2c, Ascl1, Gm10561, Fxn, Ppp1r14b, Tuba1b, Mt3
Thrsp, Apoe, Hmgn2, Hspe1, Crym, Mdk, Slc38a1, Sfrp1, Hes6, Prdx6
Ramp1, Ptgds, Txnip, Ubb, Nme2, Gpr37l1, Mt2, Gstm5, Cst3, Cmtm5
Computing nearest neighbors
Only one graph name supplied, storing nearest-neighbor graph only
Head_batch3 <- Head %>% subset(cells = batch3_lw_cells) %>%
SCTransform() %>%
RunPCA(assay = "SCT",
npcs = 100) %>%
FindNeighbors(dims = 1:10,
k.param = 30,
return.neighbor = T)
Running SCTransform on assay: RNA
Running SCTransform on layer: counts.3
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 16183 by 1238
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 1238 cells
Found 45 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 16183 genes
Computing corrected count matrix for 16183 genes
Calculating gene attributes
Wall clock passed: Time difference of 15.27602 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|==================== | 25%
|
|======================================== | 50%
|
|============================================================ | 75%
|
|================================================================================| 100%
Warning: Different cells and/or features from existing assay SCTSet default assay to SCT
PC_ 1
Positive: Cxcl10, Crym, Klf2, S100a6, Isg15, Atf3, Serpina3n, Fth1, Cebpd, H2-K1
Junb, Ccl2, H2-D1, Ier3, Ifit3, Ptx3, Cp, Sdc4, Ifitm3, Bst2
Rarres2, Vim, Mt1, Tspo, Nfkbia, Gadd45g, Mt2, Isg20, Tnfaip2, Ctsb
Negative: Nrg1, Nrxn1, Lsamp, Gm29260, Ptprz1, Dgkb, Mir99ahg, Zbtb20, Gpc5, Ntm
Pcdh9, Auts2, Msi2, Egfr, Fars2, Dcc, Ctnna2, Sparcl1, Eda, Lama2
Rora, Nlgn1, Negr1, Zfpm2, Celf2, Lrrc4c, Nol4, Cadm2, Pbx1, Dtna
PC_ 2
Positive: Cxcl10, Rsad2, Lrrc7, Ifi27l2a, Egfr, Rnf213, Marcks, Adgrl3, Iigp1, Crym
Ccnd1, Emp1, Cacnb2, Ccl2, S100a10, Slc4a7, Ptx3, Sulf1, Erbb4, Ncam2
Ext1, Npas3, Tmsb4x, Chl1, Cadps, Gbp5, Gbp2, Atrnl1, Trib1, Arid5b
Negative: Aldoc, Clu, Nrg1, Mt1, Cst3, Thbs4, Cldn10, Cmss1, Ntrk2, Junb
Gfap, Gm42418, Apoe, Gstm5, Igfbp5, Rmst, Mt3, Acsl3, Plpp3, Pla2g7
Fbxo2, Dbi, Fxyd1, Btg2, Cpe, Gadd45g, Glul, Slc1a2, Dusp1, Fam107a
PC_ 3
Positive: Cxcl10, Ifit2, Rsad2, Herc6, Gbp5, Gbp6, Ifit3, Slc1a2, Crym, Gbp2
Meg3, Nkain2, Ptprt, Slc1a3, Lsamp, Clu, Iigp1, Crot, Nrxn1, Cp
Igtp, Ifih1, Nlgn1, Trpm3, Gbp7, Pcdh9, Parp14, Serping1, Fndc3a, Isg15
Negative: Gm42418, AY036118, Cmss1, Rpl13, Rps24, Rpl32, Tmsb4x, Rpsa, Gm10260, Rplp0
Rpl23, Rpl35, Fabp7, Rps2, Rps8, Rpl15, Rps18, Rpl41, Rps12, Rps5
Rplp1, Eef1a1, Cdk8, Hist1h4d, Ppp1r14b, Hist1h1e, Marcks, Ckb, S100a10, Rpl29
PC_ 4
Positive: Mt1, Sparcl1, Gm42418, Mt2, Nrg1, Pla2g7, Cxcl10, Cmss1, Aqp4, Gbp5
Htra1, Ccl2, AY036118, Aldoc, Clu, Ifit2, Tnfaip2, Camk1d, Gfap, Psmb10
Mdga2, Gm29260, Cdk8, Cp, Hopx, Peak1, Cxcl14, Kctd1, Tspo, Id2
Negative: Klf2, Btg2, Hes5, Dnajb1, Klf4, Ier2, Atf3, Junb, Jun, Egr1
Fos, Egr2, Homer1, Dusp1, Zfp36, Ddit3, Fosb, Egr3, Nfatc2, Ccn1
Klf6, Hes1, Rhob, Meg3, Jund, Creb5, Vcl, Mlf1, Nr4a1, Gadd45g
PC_ 5
Positive: Hes5, Gm42418, Isg15, Rsad2, Iigp1, Fabp7, Cmss1, Ifit2, Pdzph1, Ifit3
AY036118, Cxcl10, Crym, Gm4951, Ifi27l2a, Ifit1, AW011738, Ifit3b, Ifi44, Ifih1
Hes1, Cmpk2, Ckb, Egfr, Parp11, Igtp, Ndufc2, Ier2, Rnf213, Cdk8
Negative: Ptx3, Kcnip4, Nrg1, Emp1, Trib1, Mt1, Csmd1, S100a10, Mt2, Ncam2
Dlg2, Vim, Gfap, Anxa2, Cadm2, Serpina3n, Hspb1, Clu, Sorcs1, Sdc4
Mrps6, Chl1, Gap43, Ifrd1, Dhrs1, S100a11, Msn, Pdpn, Gm4876, Atf3
Computing nearest neighbors
Only one graph name supplied, storing nearest-neighbor graph only
#Inferring unlabeled b1 cells in batch 2:
wpre_index_batch2 <- which(Head_batch2@neighbors$SCT.nn@cell.names %in% wpre_cells_batch2)
good.neighbors_batch2 <- c()
bad.neighbors <- wpre_index_batch2 #Excluding labeled cells from the list of neighbors
for (i in wpre_index_batch2) {
all.neighbors <- Head_batch2@neighbors$SCT.nn@nn.idx[i,2:30]
good.neighbors_batch2 <- c(good.neighbors_batch2,
(all.neighbors[!all.neighbors %in% bad.neighbors])[1:3])
bad.neighbors <- c(bad.neighbors, good.neighbors_batch2) %>% unique()
}
wpre_neighbors_batch2 <- Head_batch2@neighbors$SCT.nn@cell.names[good.neighbors_batch2]
wpre_neighbors_batch2 <- wpre_neighbors_batch2[!is.na(wpre_neighbors_batch2)]
#Inferring unlabeled b1 cells in batch 3:
wpre_index_batch3 <- which(Head_batch3@neighbors$SCT.nn@cell.names %in% wpre_cells_batch3)
good.neighbors_batch3 <- c()
bad.neighbors <- wpre_index_batch3
for (i in wpre_index_batch3) {
all.neighbors <- Head_batch3@neighbors$SCT.nn@nn.idx[i,2:30]
good.neighbors_batch3 <- c(good.neighbors_batch3, (setdiff(all.neighbors, bad.neighbors))[1:3])
bad.neighbors <- c(bad.neighbors, good.neighbors_batch3)
}
wpre_neighbors_batch3 <- Head_batch3@neighbors$SCT.nn@cell.names[good.neighbors_batch3]
wpre_neighbors_batch3 <- wpre_neighbors_batch3[!is.na(wpre_neighbors_batch3)]
wpre.data[, "tdtom"] <- "TdTomato-"
wpre.data[c(wpre_neighbors_batch2, wpre_neighbors_batch3),"tdtom"] <- "TdTomato+ NN"
wpre.data[c(wpre_cells_batch2, wpre_cells_batch3), "tdtom"] <- "TdTomato+"
wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch, tdtom) %>%
summarise(n = n()) %>%
left_join(wpre.data %>%
filter(region != "Wedge") %>%
group_by(batch) %>%
summarize(batch_n = n()), by = "batch") %>%
mutate(pct = n/batch_n)
`summarise()` has grouped output by 'batch'. You can override using the `.groups` argument.
wpre.data %>%
filter(region != "Wedge" & batch == 'batch_2') %>%
group_by(age, tdtom) %>%
summarise(n = n()) %>%
left_join(wpre.data %>%
filter(region != "Wedge" & batch == 'batch_2') %>%
group_by(age) %>%
summarize(age_n = n(), by = "age")) %>%
mutate(pct = n/age_n)
`summarise()` has grouped output by 'age'. You can override using the `.groups` argument.Joining with `by = join_by(age)`
Head$tdtom <- wpre.data$tdtom
Head@meta.data %>%
rownames_to_column('bc') %>%
left_join(wpre.data %>% rownames_to_column('bc') %>% select(bc, tdtom), by = 'bc') %>%
column_to_rownames('bc')
Head@meta.data[, "b_type"] <- "b2"
Head@meta.data[c(wpre_cells_batch2, wpre_cells_batch3, wpre_neighbors_batch2, wpre_neighbors_batch3),"b_type"] <- "b1"
mDimPlot(Head, group.by = "b_type", shuffle = T)
mDimPlot(Head, group.by = "b_type", split.by = "age")
Head@active.assay = "RNA"
non_infected_features <- rownames(Head)[!rownames(Head) %in% c('Wpre', 'GFP', 'Cre')]
Head_batch1 <- subset(Head, subset = batch == "batch_1", features = non_infected_features) %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA(assay = "RNA",
npcs = 20) %>%
FindNeighbors(dims = 1:20)
Normalizing layer: counts.1
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.1
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Centering and scaling data matrix
|
| | 0%
|
|======================================== | 50%
|
|================================================================================| 100%
PC_ 1
Positive: Egfr, Kcnh7, Sox11, Lrrfip1, Slc4a7, Tcf12, Tmsb4x, Map2, Setbp1, Ppp1r14b
Hnrnpab, Nrxn3, Sipa1l1, Shc3, Dpysl3, Ccnd2, Rpl32, Sv2c, Magi1, Dpp6
Chd7, Marcks, Nell2, Nfib, Itgb1, Hnrnpa1, Fgd4, Slit2, Auts2, Ybx1
Negative: Clu, Mt1, Mt3, Fxyd1, Mt2, S100a6, Rgcc, Junb, Ccn1, Sparc
Id2, Pla2g7, Id3, Cebpb, Gfap, Thbs4, Rarres2, Sat1, Gsn, Nap1l5
Serpine2, Zfp36, Bhlhe40, Gm39325, Gem, Dusp1, Ramp1, Dad1, Gadd45g, Slc25a3
PC_ 2
Positive: Riiad1, Crym, Fam183b, Basp1, Cebpd, Dpyd, Ncam2, Ramp1, Pifo, 1700094D03Rik
Klf2, Cfap126, Fas, Capsl, Foxj1, Adgrl3, Cfap43, Rsph1, Ifitm3, C4b
Dnah9, S100a6, Slc38a1, Dnah6, Mns1, Nnat, Mlf1, C230072F16Rik, Meig1, Sorbs1
Negative: Nrg1, Urah, Gpc5, Cdk11b, Dclk1, Gm39325, Vps37b, Dio2, Pla2g7, Thrsp
Slc25a3, Epha5, Snhg15, Trpm3, Gm26512, Gsx2, Gm29260, Ppp1r3g, Ptprz1, Sgip1
Nkain2, Usp2, Emd, Dcc, Gm10563, Pls1, Pde10a, Arih1, Mpp7, Pmepa1
PC_ 3
Positive: Stmn2, Dlx6os1, Tubb3, Gm38505, Igfbpl1, Gad2, Dcx, Nav3, Sp9, Celf4
Adarb2, Cacna1c, Myt1l, Sp8, Dlx5, Abracl, 2610307P16Rik, Myt1, Edil3, Anks1b
Shtn1, Bcl11b, Nsg2, Brinp2, Dlx1, Mir124a-1hg, Gad1, Crmp1, Elavl4, Arx
Negative: Egr3, Notum, Fabp7, Slit2, Egfr, Creb5, Nfatc2, Nptx2, Rfx4, Rasal2
Slco1c1, Ddit3, Hspa5, Hmgcs1, Adgrl3, Sik2, Aspscr1, Prss56, Crym, Slc1a3
Shc3, Egr2, Emp1, Pde10a, Dnajc1, Zfp521, Ascl1, Dnajb1, Ankrd28, Dusp10
PC_ 4
Positive: Maff, Nr4a1, Slc25a3, Arih1, Cebpb, Nr4a3, Nop53, Vps37b, Usp2, Emd
Lmna, Pde10a, Btbd9, Braf, Homer1, Ifrd1, Atf3, Tpm4, Dad1, Nr4a2
Cdk11b, Zfp36, Anxa2, Dclk1, Srxn1, Ptma, Ppp1r15a, Isy1, Hells, Ccn1
Negative: Grm3, Ntm, Pcdh9, Cntnap2, Ptprt, Sox2ot, Rmst, Gm4258, B3galt1, Trim9
Dcc, Fjx1, Frmd4a, Veph1, Gm29260, Gm35552, Ptprz1, Hdac9, Auts2, Lrrc4c
Grid2, Lama2, Kctd1, 9330159F19Rik, Fabp7, Tfap2c, Pter, Zfpm2, Sall3, Cdh10
PC_ 5
Positive: Arih1, Braf, Usp2, Zswim6, Meis2, Vps37b, Slco1c1, Nr4a1, Ctnna2, Btbd9
Maff, Pak1, Pde10a, Cdk11b, Slc1a3, Nr4a3, Homer1, Dclk1, Camk2d, Ccn1
Zdbf2, Tubb2a, Pcdh9, Nr4a2, Atxn7, Creb5, Meg3, Emd, Clasp2, Sik2
Negative: Hells, Chaf1b, Mcm6, Mcm5, Ung, Uhrf1, Pclaf, Mcm3, Cenpm, Dut
Mcm2, Cenph, Gins2, E2f1, Tcf19, Lig1, Gmnn, Clspn, Gadd45g, Dctpp1
Dhfr, Dtl, Nme2, Fignl1, Rmi2, Anp32b, Mcm7, Spc24, H2afz, Top2a
Warning: Number of dimensions changing from 100 to 20Computing nearest neighbor graph
Computing SNN
Head_batch23 <- subset(Head, subset = batch %in% c("batch_2", "batch_3"))
Head_batch23 <- Head_batch23 %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA(assay = "RNA",
npcs = 20)
Normalizing layer: counts.2
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.3
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.2
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.3
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Centering and scaling data matrix
|
| | 0%
|
|======================================== | 50%
|
|================================================================================| 100%
PC_ 1
Positive: Bst2, Tspo, Fth1, Ifit3, Psmb8, H2-K1, H2-D1, Isg15, Psmb10, Ifit3b
Serping1, Clu, Ifitm3, Gbp2, Cxcl10, S100a6, B2m, H2-Q4, Cebpd, Serpina3n
Gfap, Sdc4, H2-T23, Rarres2, Isg20, C1ra, Gm42418, Irf7, Gbp5, Ly6e
Negative: Auts2, Dgkb, Pbx1, Nrxn1, Sox11, Ptprz1, Kcnh7, Nol4, Lsamp, Egfr
Magi2, Pak3, Map2, Eda, Slit2, Zeb1, Marcks, Ccnd2, Dscaml1, Ctnna2
Baz2b, Grin2b, Ntm, Chd7, Setbp1, Nell2, Sall3, Cdh4, Maml3, Coro1c
PC_ 2
Positive: Crym, Cxcl10, Spon1, Klf6, Atf3, Rsad2, Klf2, Egr1, Ascc3, Six3
Nfatc2, Tagln2, Emp1, Zfp36, Jun, Ddit3, Egr3, Gbp2, Egr2, Samd9l
Lmna, Crot, Ext1, Dnajb1, Gbp5, Notum, Iigp1, Trib1, Irf7, Efhb
Negative: Pla2g7, Mt1, Sparcl1, Thbs4, Nrg1, Fxyd6, Clu, Gm42418, Kctd1, Gfap
Rmst, Gm29260, Rlbp1, Gsn, Gabrb1, 9530026P05Rik, Glycam1, Epha5, Hopx, Mdga2
Mt2, Slc15a2, Nr2f1, Prex2, AY036118, A330093E20Rik, Pls1, Gm35552, Garnl3, Dcc
PC_ 3
Positive: Gbp7, Anks1b, Rnf213, Iigp1, Znfx1, Ptprz1, Parp14, Oasl2, Gbp5, Atp10a
Rbms3, Ifit2, Ddx58, Trim30a, Mdga2, Gm29260, Parp11, Herc6, Slfn8, Igtp
Sparcl1, Ifih1, Fndc3a, Ifit3, Psmb10, Angpt1, Gbp6, Csmd2, Tor3a, Gm4951
Negative: Btg2, Klf4, Dusp1, Junb, Dnajb1, Zfp36, Ier2, Jun, Klf2, Mlf1
Maff, Nr4a1, Egr1, Egr2, Ddit3, Cebpb, Gm26887, Atf3, 6330403K07Rik, Homer1
Meg3, Foxj1, Ccn1, Gadd45g, Notum, Klf6, Hes1, Ubald1, Idi1, Insig1
PC_ 4
Positive: Tmsb4x, Ppp1r14b, AY036118, Tmsb10, H2afz, Gm42418, Tubb5, Hmgn2, Lrrc7, Selenoh
Lrrfip1, S100a10, Top2a, Stmn1, Hells, Cdca8, Cks2, Slit1, Uhrf1, Mki67
Cenpk, Mcm2, Ccnd1, Gm38505, Dcx, Nsg2, Gm10260, Mcm5, Cenpf, Dpysl3
Negative: Lsamp, Nrxn1, Nkain2, Cadm2, Csmd1, Gabrb1, Pcdh9, Nlgn1, Grid2, Prex2
Rgs6, Ppp2r2b, Nrg1, Herc6, Dlgap1, Luzp2, Dgkb, Slc35f1, Vcl, Bmpr1b
Dcc, Nckap5, Timp3, Hdac8, Auts2, Klf4, Ntm, Lrrc4c, Plagl1, Fgf13
PC_ 5
Positive: Fabp7, Egfr, Lrrfip1, Ptprz1, Olig2, Nptx2, Hes5, Sv2c, Marcks, Sall3
Thrsp, Ascl1, Slc4a7, Ccnd1, Rhcg, Shc3, Olig1, Adam12, Sema5b, Cdk6
Dll1, Crlf1, Me3, A930037H05Rik, Ifi27l2a, Nes, Pvt1, Hs2st1, Rnf213, Coro1c
Negative: Dlx6os1, Dcx, Nav3, Gm38505, Stmn2, Tubb3, Cacna1c, Sp9, Clu, Shtn1
Gad2, Adarb2, Myt1l, Celf4, Igfbpl1, Islr2, Gad1, Pak7, Synpr, Kcnh8
Mt1, Tenm2, Pcdha-all, Dlx5, Bcl11b, Plxna4, Enox1, Basp1, Nrxn3, Stmn4
Warning: Number of dimensions changing from 100 to 20
Head_batch23 <- IntegrateLayers(Head_batch23, method = CCAIntegration, orig.reduction = "pca", dims = 1:20)
Finding all pairwise anchors
| | 0 % ~calculating
Running CCA
Merging objects
Finding neighborhoods
Finding anchors
Found 2139 anchors
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s
Merging dataset 1 into 2
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data
Warning: Number of dimensions changing from 100 to 20
Head.anchors <- Head_batch23 %>%
FindTransferAnchors(reference = .,
query = Head_batch1,
dims = 1:20,
reference.reduction = "integrated.dr")
Projecting cell embeddings
Finding neighborhoods
Finding anchors
Found 1308 anchors
predictions <- TransferData(anchorset = Head.anchors, refdata = Head_batch23$b_type, dims = 1:20)
Finding integration vectors
Finding integration vector weights
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels
inferredb1_batch1_cells <- predictions %>% filter(predicted.id != "b2") %>% rownames()
Head@meta.data[inferredb1_batch1_cells, "tdtom"] <- "TdTomato+ LT"
Head@meta.data[inferredb1_batch1_cells, "b_type"] <- "b1"
#---------------------------------
# Head@active.assay = "RNA"
# non_infected_features <- rownames(Head)[!rownames(Head) %in% c('Wpre', 'GFP', 'Cre')]
#
# Head_batch1 <- Head %>%
# subset(subset = batch == "batch_1", features = non_infected_features) %>%
# NormalizeData() %>%
# FindVariableFeatures() %>%
# ScaleData() %>%
# RunPCA(assay = "RNA",
# npcs = 20) %>%
# FindNeighbors(dims = 1:20)
#
#
# Head_batch23 <- Head %>% subset(subset = batch != "batch_1", features = non_infected_features)
# Head_batch23 <- JoinLayers(Head_batch23, assay = "RNA")
# Head_batch23[["RNA"]] <- split(Head_batch23[["RNA"]], f = Head_batch23$batch)
#
# Head_batch23 <- Head_batch23 %>%
# NormalizeData() %>%
# FindVariableFeatures() %>%
# ScaleData() %>%
# RunPCA(assay = "RNA",
# npcs = 20)
# Head_batch23 <- IntegrateLayers(Head_batch23, method = CCAIntegration, orig.reduction = "pca", dims = 1:20)
#
# Head.anchors <- Head_batch23 %>%
# FindTransferAnchors(reference = .,
# query = Head_batch1,
# dims = 1:20,
# reference.reduction = "pca")
#
# predictions <- TransferData(anchorset = Head.anchors,
# refdata = Head_batch23$b_type,
# dims = 1:20)
#
# predictions %>% filter(predicted.id =='b1')
#
# inferredb1_batch1_cells <- predictions %>% filter(predicted.id == "b1") %>% rownames()
# Head@meta.data[inferredb1_batch1_cells, "tdtom"] <- "TdTomato+ LT"
# Head@meta.data[inferredb1_batch1_cells, "b_type"] <- "b1"
Head@active.assay <- 'SCT'
mFeaturePlot(Head, features = c("S100a6", "Egfr", "Crym", "Urah", "Tfap2c"), order = T, ncol = 5, legend = F)
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
mDimPlot(Head, group.by = "age", order = T) +
scale_color_manual(values = c('grey90', 'tomato'))
mDimPlot(Head, group.by = "activation_state", split.by = "age", shuffle = T) +
scale_color_manual(values = c('#1c9099', '#a6bddb'), name = 'Activation State')
mFeaturePlot(Head, features = c("Crym", "Egfr"), split.by = "age", order = T)
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
mDimPlot(Head, group.by = "b_type", shuffle = T) +
scale_color_manual(values = c("b1" = "skyblue1", "b2" = "mediumblue"), labels = c('B1', 'B2'))
mDimPlot(Head, group.by = "b_type", split.by = "age", shuffle = T) +
scale_color_manual(values = c("b1" = "skyblue1", "b2" = "mediumblue"), labels = c('B1', 'B2'))
mDimPlot(Head, group.by = "tdtom", order = T) +
scale_color_manual(values = c("TdTomato+" = "magenta",
"TdTomato+ NN" = "darkorange",
"TdTomato+ LT" = 'darkgreen'),
na.value = 'grey90')
mDimPlot(Head, group.by = "tdtom", split.by = 'tdtom', order = T) +
scale_color_manual(values = c("TdTomato+" = "magenta",
"TdTomato+ NN" = "darkorange",
"TdTomato+ LT" = 'darkgreen'),
na.value = 'grey90')
mDimPlot(Head, group.by = "SCT_snn_res.1", shuffle = T, label = T, legend = F)
mDimPlot(Head, group.by = "SCT_snn_res.1", split.by = 'region', shuffle = T, label = T, legend = F)
Head@meta.data %>%
ggplot(aes(x = age, fill = b_type)) +
geom_bar(position = 'fill') +
theme_classic() +
scale_fill_manual(values = c("b1" = "skyblue1", "b2" = "mediumblue"), labels = c('B1', 'B2')) +
scale_y_continuous(labels = scales::percent, expand = c(0,0)) +
labs(x = "Age", y = "Fraction of cells", fill = "B cell type")
ggsave('../figures/b1b2_barplot.pdf', width = 5, height = 3)
a <- Head@meta.data %>%
filter(region == 'LW' & age == 'p30') %>%
ggplot(aes(x = b_type, fill = activation_state)) +
geom_bar() +
theme_classic() +
scale_y_continuous(expand = c(0, 0)) +
scale_x_discrete(labels = c('B1', 'B2')) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5)) +
scale_fill_manual(values = c('#1c9099', '#a6bddb'), name = 'Activation State') +
labs(x = NULL, y = 'Count')
b <- Head@meta.data %>%
filter(region == 'LW' & age == 'p30') %>%
ggplot(aes(x = b_type, fill = activation_state)) +
geom_bar(position = 'fill') +
theme_classic() +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
scale_x_discrete(labels = c('B1', 'B2')) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5)) +
scale_fill_manual(values = c('#1c9099', '#a6bddb'), name = 'Activation State') +
labs(x = NULL, y = 'Share')
a + b + plot_layout(guides = 'collect') + plot_annotation(caption = 'Only cells in the LW region at p30 are shown.')
c <- Head@meta.data %>%
filter(region == 'LW') %>%
ggplot(aes(x = age, fill = activation_state)) +
geom_bar(position = 'fill') +
theme_classic() +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5)) +
scale_fill_manual(values = c('#1c9099', '#a6bddb'), name = 'Activation State') +
labs(x = NULL, y = 'Share')
b + c + plot_layout(guides = 'collect', axis_titles = 'collect_y') + plot_annotation(caption = 'Only cells in the LW region are shown.')
ggsave("../figures/activation_states.pdf", width = 5, height = 3)
#- DEGs
Idents(Head) <- "b_type"
HeadQ <- Head %>% subset(subset = activation_state == 'quiescent') %>%
SCTransform(return.only.var.genes = F) %>%
PrepSCTFindMarkers()
Running SCTransform on assay: RNA
Running SCTransform on layer: counts.1
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 17603 by 2321
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 2321 cells
Found 100 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 17603 genes
Computing corrected count matrix for 17603 genes
Calculating gene attributes
Wall clock passed: Time difference of 26.05989 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|==== | 4%
|
|======== | 8%
|
|============= | 12%
|
|================= | 17%
|
|===================== | 21%
|
|========================= | 25%
|
|============================= | 29%
|
|================================== | 33%
|
|====================================== | 38%
|
|========================================== | 42%
|
|============================================== | 46%
|
|================================================== | 50%
|
|======================================================= | 54%
|
|=========================================================== | 58%
|
|=============================================================== | 62%
|
|=================================================================== | 67%
|
|======================================================================== | 71%
|
|============================================================================ | 75%
|
|================================================================================ | 79%
|
|==================================================================================== | 83%
|
|======================================================================================== | 88%
|
|============================================================================================= | 92%
|
|================================================================================================= | 96%
|
|=====================================================================================================| 100%
Running SCTransform on layer: counts.2
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 14890 by 323
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 323 cells
Found 42 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 14890 genes
Computing corrected count matrix for 14890 genes
Calculating gene attributes
Wall clock passed: Time difference of 6.948273 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|===== | 5%
|
|========== | 10%
|
|=============== | 15%
|
|==================== | 20%
|
|========================= | 25%
|
|============================== | 30%
|
|=================================== | 35%
|
|======================================== | 40%
|
|============================================= | 45%
|
|================================================== | 50%
|
|======================================================== | 55%
|
|============================================================= | 60%
|
|================================================================== | 65%
|
|======================================================================= | 70%
|
|============================================================================ | 75%
|
|================================================================================= | 80%
|
|====================================================================================== | 85%
|
|=========================================================================================== | 90%
|
|================================================================================================ | 95%
|
|=====================================================================================================| 100%
Running SCTransform on layer: counts.3
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 16617 by 1889
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 1889 cells
Found 48 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 16617 genes
Computing corrected count matrix for 16617 genes
Calculating gene attributes
Wall clock passed: Time difference of 16.08269 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|==== | 4%
|
|========= | 9%
|
|============= | 13%
|
|================== | 17%
|
|====================== | 22%
|
|========================== | 26%
|
|=============================== | 30%
|
|=================================== | 35%
|
|======================================== | 39%
|
|============================================ | 43%
|
|================================================ | 48%
|
|===================================================== | 52%
|
|========================================================= | 57%
|
|============================================================= | 61%
|
|================================================================== | 65%
|
|====================================================================== | 70%
|
|=========================================================================== | 74%
|
|=============================================================================== | 78%
|
|=================================================================================== | 83%
|
|======================================================================================== | 87%
|
|============================================================================================ | 91%
|
|================================================================================================= | 96%
|
|=====================================================================================================| 100%
Centering data matrix
|
| | 0%
|
|==== | 4%
|
|======== | 8%
|
|============= | 12%
|
|================= | 17%
|
|===================== | 21%
|
|========================= | 25%
|
|============================= | 29%
|
|================================== | 33%
|
|====================================== | 38%
|
|========================================== | 42%
|
|============================================== | 46%
|
|================================================== | 50%
|
|======================================================= | 54%
|
|=========================================================== | 58%
|
|=============================================================== | 62%
|
|=================================================================== | 67%
|
|======================================================================== | 71%
|
|============================================================================ | 75%
|
|================================================================================ | 79%
|
|==================================================================================== | 83%
|
|======================================================================================== | 88%
|
|============================================================================================= | 92%
|
|================================================================================================= | 96%
|
|=====================================================================================================| 100%
Centering data matrix
|
| | 0%
|
|===== | 5%
|
|========== | 10%
|
|=============== | 15%
|
|==================== | 20%
|
|========================= | 25%
|
|============================== | 30%
|
|=================================== | 35%
|
|======================================== | 40%
|
|============================================= | 45%
|
|================================================== | 50%
|
|======================================================== | 55%
|
|============================================================= | 60%
|
|================================================================== | 65%
|
|======================================================================= | 70%
|
|============================================================================ | 75%
|
|================================================================================= | 80%
|
|====================================================================================== | 85%
|
|=========================================================================================== | 90%
|
|================================================================================================ | 95%
|
|=====================================================================================================| 100%
Centering data matrix
|
| | 0%
|
|==== | 4%
|
|========= | 9%
|
|============= | 13%
|
|================== | 17%
|
|====================== | 22%
|
|========================== | 26%
|
|=============================== | 30%
|
|=================================== | 35%
|
|======================================== | 39%
|
|============================================ | 43%
|
|================================================ | 48%
|
|===================================================== | 52%
|
|========================================================= | 57%
|
|============================================================= | 61%
|
|================================================================== | 65%
|
|====================================================================== | 70%
|
|=========================================================================== | 74%
|
|=============================================================================== | 78%
|
|=================================================================================== | 83%
|
|======================================================================================== | 87%
|
|============================================================================================ | 91%
|
|================================================================================================= | 96%
|
|=====================================================================================================| 100%
Warning: Different cells and/or features from existing assay SCTSet default assay to SCT
Found 3 SCT models. Recorrecting SCT counts using minimum median counts: 6415
| | 0 % ~calculating
|+++++++++++++++++ | 33% ~38s
|++++++++++++++++++++++++++++++++++ | 67% ~12s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=40s
btype_markers <- FindAllMarkers(HeadQ,
logfc.threshold = 0,
min.pct = 0,
only.pos = T)
Calculating cluster b1
Calculating cluster b2
btype_markers <- btype_markers %>%
group_by(cluster) %>%
mutate(rank = rank(p_val_adj, ties.method = "first")) %>%
arrange(cluster, rank)
#- Volcano Plot
validated.genes = c("Atf3", "Ptprz1", "Riiad1", "FoxJ1", "Gadd45b", "Zeb1", "Tagln2", "Emp1", "Anxa2")
regional.genes = c("Crym", "Nrg1", "Klf2", "Cebpd", "Gm29260", "Pax6", "Rlbp1", "Nkx6-2")
b1_degenes = btype_markers %>% filter(cluster == "b1") %>% rownames()
b2_degenes = btype_markers %>% filter(cluster == "b2") %>% rownames()
genes_to_highlight = c(b1_degenes[1:5], b2_degenes[1:5],validated.genes, regional.genes)
volcano.data = btype_markers %>%
mutate(labs = if_else(gene %in% genes_to_highlight, gene, "")) %>%
mutate(significant = case_when(p_val_adj < 0.05 &
cluster == "b1" ~ "b1",
p_val_adj < 0.05 & cluster == "b2" ~ "b2",
TRUE ~ "Not differentially expressed")) %>%
mutate(significant = if_else(gene%in%regional.genes, "regional.genes", significant)) %>%
mutate(log2FC = case_when(cluster == "b2" ~ avg_log2FC*-1,
TRUE ~ avg_log2FC)) %>%
mutate(significant = if_else(gene%in%validated.genes, "Validated gene", significant)) %>%
mutate(priority = gene %in% genes_to_highlight) %>%
arrange(priority, gene) %>%
select(-priority)
volcano.data %>%
ggplot(aes(log2FC, -log(base = 10, p_val_adj), label = labs)) +
geom_point(aes(col = significant), stroke = 0, size = 4, alpha = 0.5) +
scale_color_manual(name = NULL, labels = c("Upregulated",
"Downregulated",
"Not differentially expressed",
"Regionalization-associated genes",
"Validated gene"
), values = c("skyblue1", "mediumblue", "grey75","darkgreen","magenta")) +
geom_text_repel(max.overlaps = Inf, size=5, box.padding=0.3) +
theme_classic() +
theme(panel.grid.minor = element_blank(),
text = element_text(family = "Helvetica", size = 15)
) +
labs(title = "Differential expression between B1 and B2 cells",
x = "Average log2FC",
y = "-log10 adjusted P value")
ggsave("data/btype_volcano.pdf", width = 6, height = 4)
#-Cell Cycle Score
g2m.genes
[1] "Hmgb2" "Cdk1" "Nusap1" "Ube2c" "Birc5" "Tpx2" "Top2a" "Ndc80" "Cks2"
[10] "Nuf2" "Cks1b" "Mki67" "Tmpo" "Cenpf" "Tacc3" "Fam64a" "Smc4" "Ccnb2"
[19] "Ckap2l" "Ckap2" "Aurkb" "Bub1" "Kif11" "Anp32e" "Tubb4b" "Gtse1" "Kif20b"
[28] "Hjurp" "Cdca3" "Hn1" "Cdc20" "Ttk" "Cdc25c" "Kif2c" "Rangap1" "Ncapd2"
[37] "Dlgap5" "Cdca2" "Cdca8" "Ect2" "Kif23" "Hmmr" "Aurka" "Psrc1" "Anln"
[46] "Lbr" "Ckap5" "Cenpe" "Ctcf" "Nek2" "G2e3" "Gas2l3" "Cbx5" "Cenpa"
ggsave("..figures/activation_states.pdf", width = 5, height = 3)
Cannot find directory ]8;;file:///media/data5/marcos/b2/manuscript/2.analysis/..figures..figures]8;;.
ℹ Would you like to create a new directory?
1: Yes
2: No
2
Error in `ggsave()`:
! Cannot find directory ]8;;file:///media/data5/marcos/b2/manuscript/2.analysis/..figures..figures]8;;.
ℹ Please supply an existing directory or use `create.dir = TRUE`.
Backtrace:
1. ggplot2::ggsave("..figures/activation_states.pdf", width = 5, height = 3)
Head@active.assay <- "SCT"
b.markers.list = list(b1 = c("Atf3", "Riiad1", "Foxj1", "Gadd45b", "Tagln2", "Emp1"))
Head = AddModuleScore(Head, features = b.markers.list)
Head@meta.data$B1_score <- NULL
Head@meta.data = Head@meta.data %>% dplyr::rename(B1_score = Cluster1)
#Re-scaling Module Scores:These lines rescale the module scores (B1_score and B2_score) so that they have a mean of 0 and a standard deviation of 1. This step standardizes the scores, making them comparable across different modules.
Head$B1_score = Head$B1_score %>% rescale()
mFeaturePlot(Head, features = "B1_score", order = T) &
scale_color_distiller(type = 'seq', palette = 16, direction = 1)
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
saveRDS(Head, "data/Head_withp365.rds")
mFeaturePlot(Head, features = c('Atf3', 'Tagln2', 'Emp1'), order = T, ncol = 3, legend = F)
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
mFeaturePlot(HeadQ,
features = b.markers.list$b1,
split.by = 'b_type',
order = T,
legend = F)
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
DotPlot(HeadQ, b.markers.list$b1, group.by = 'b_type', scale = F )
Head@meta.data %>% filter(activation_state == 'quiescent') %>%
ggplot(aes(x = age, y = B1_score, color = b_type)) +
geom_jitter(height = 0, alpha = 0.2, size = 0.5) +
geom_boxplot(outliers = F, fill = NA, size = 0.8) +
theme_classic() +
scale_y_continuous() +
scale_color_manual(name = 'B Type',
values = c('b1' = 'skyblue', 'b2' = 'mediumblue'),
labels = c('B1', 'B2')) +
labs(title = 'B1 Score in quiescent B cells', x = 'Age', y = 'B1 Score') +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
ggsave('../figures/b1_scores_by_age.pdf')
Saving 7.29 x 4.51 in image
#-Heatmap
#Main color
col_values <- seq(-2, 2, 0.5)
b1_score_values <- seq(0, 1, 0.125)
col_fun <- colorRamp2(c(-2, 0, 2), c("cyan", "white", "magenta"))
#annotations colors
top_anno_colors <- list(Dissection = c("LV" = "grey90", "LW" = "purple", "Wedge" = "forestgreen"),
`B Cell type` = c("b1" = 'skyblue1', "b2" = "mediumblue"),
`B1 Score` = colorRamp2(b1_score_values, brewer.pal(n = length(b1_score_values), name = "YlGnBu")))
pdf("../figures/b1b2_degenes_heatmap.pdf", width = 8.2, height = 5.8)
hm
dev.off()
null device
1
#P9+P30+P365 ——————————-
# Head <- subset(integrated_exp, cells = b_cells)
Head <- Head %>%
SCTransform() %>%
RunPCA(assay = "SCT", npcs = 100) %>%
IntegrateLayers(method = CCAIntegration,
normalization.method = "SCT",
verbose = F) %>%
RunUMAP(dims = 1:100, reduction = "integrated.dr") %>%
FindNeighbors(reduction = "integrated.dr", dims = 1:100) %>%
FindClusters(resolution = seq(0.5, 2, 0.5), graph.name = "SCT_snn")
DimPlot(Head, label=T, shuffle = T, group.by ="SCT_snn_res.1") +
NoLegend() +
coord_fixed() +
NoAxes()